%==========================================================================
%=  This file is part of the Sound Restoration Project
%=  (c) Copyright Industrial Mathematics Institute
%=                University of South Carolina, Department of Mathematics
%=  ALL RIGHTS RESERVED
%=
%=  Author: Borislav Karaivanov
%==========================================================================

%==========================================================================
% List of the files on which this procedure depends:
%
% findBestFitToSignal.m
% fitReduceFitSignal.m
% visualizeSoundCalibrationSingleMask.m
%
% curve fitting toolbox
%
%==========================================================================

%==========================================================================
% The function "calibrateSignalsSingleMask" takes an array of signals of
% the same length, possibly reduces them in length by discarding fixed
% segments of samples at the beginning and the end of each of them, and
% calibrates them with a calibration mask based on the average of all of
% them (smoothed or fitted by a low-degree polynomial).
% INPUT: "origSignalsCellArr" is a cell column vector with cells that are
% 2D arrays whose columns hold the individual sound signals.
% "sweeps" is an optional non-negative integer specifying the number of
% smoothing passes to be applied to the raw averages of the signals.
% "span" is an optional, positive, odd integer specifying the span of the
% stencil to be used for moving average smoothing of the average signal.
% Value of 1 results in no smoothing at all.
% "degOfPolyFitArr" is an optional array of non-negative integers
% specifying the degrees of the polynomials to be fit to the average
% signal. The best fitting (in uniform metric sense) of these polynomial is
% to be sampled and used as a calibration mask. If this parameter is
% specified (and non-empty), then the values of the previous two parameters
% are ignored.
% "samplesToIgnore" is an optional, non-negative integer specifying the
% maximal allowed cumulative number of samples that can be removed at the
% beginning or at the end of the given signal in order to produce a reduced
% signal with a better polynomial fit.
% OUTPUT: "calibSignalsArr" returns a 2D array of the same width as the
% input array of signals and not greater height. Its columns hold the
% individual, (possibly reduced), calibrated signals.
% "maskArr" returns a column vector of length equal to the height of the
% array of calibrated signals. It holds the calibration mask by which the
% reduced original signals are to be divided in order to produce their
% calibrated versions.
% "leftSamplesToIgnore" returns the number of samples removed at the left
% end of the given signals.
% "rightSamplesToIgnore" returns the number of samples removed at the right
% end of the given signals.
%==========================================================================
function [calibSignalsArr, maskArr, leftSamplesToIgnore, rightSamplesToIgnore] = ...
    calibrateSignalsSingleMask(origSignalsCellArr, sweeps, span, ...
    degOfPolyFitArr, samplesToIgnore)

if (nargin < 2)
    sweeps = 1;
end
if (nargin < 3)
    span = 15;
end
if (nargin < 4)
    degOfPolyFitArr = (1:4);
end
if (nargin < 5)
    samplesToIgnore = 100;
end

    
%=========================================================================%
%========== Here is the essential portion of this function. ==============%
%=========================================================================%
% Get the common length of all signals.
signalLength = size(origSignalsCellArr{1}, 1);
% Get the number of signals.
numSignals = sum(cellfun(@(x) size(x, 2), origSignalsCellArr));

% Compute the average of the signal.
aveArr = sum(cell2mat((cellfun(@(x) mean(x, 2)*(size(x, 2)/numSignals), ...
    origSignalsCellArr, 'UniformOutput', false)')), 2);

% Find a good, smoother fit to the average signal.
[maskArr, leftSamplesToIgnore, rightSamplesToIgnore, degOfPolyFit] = ...
    findBestFitToSignal(aveArr, sweeps, span, degOfPolyFitArr, samplesToIgnore);

% Allocate memory for the calibrated signals.
calibSignalsArr = zeros(signalLength - leftSamplesToIgnore - rightSamplesToIgnore, numSignals);

% Calibrate the signals using the computed calibration mask.
ind = 0;
for cellInd = 1:numel(origSignalsCellArr)
    for columnInd = 1:size(origSignalsCellArr{cellInd}, 2)
        ind = ind + 1;
        calibSignalsArr(:, ind) = origSignalsCellArr{cellInd} ...
            ((leftSamplesToIgnore + 1):(end - rightSamplesToIgnore), columnInd)./maskArr;
    end
end
%=========================================================================%

% % Fix a batch from which a signal is to be drawn.
% batchInd = randi(numel(origSignalsCellArr));
% % Fix a signal to be visualized.
% signalToShowInd = randi(size(origSignalsCellArr{batchInd}, 2));
% % Visualize the full-length signals and save snapshots, if desired.
% visualizeSoundCalibrationSingleMask(origSignalsCellArr{batchInd}(:, signalToShowInd), ...
%     aveArr, sweeps, span, degOfPolyFit);
% % Visualize the reduced signals and save snapshots, if desired.
% visualizeSoundCalibrationSingleMask(origSignalsCellArr{batchInd}...
%     ((leftSamplesToIgnore + 1):(end - rightSamplesToIgnore), signalToShowInd), ...
%     aveArr((leftSamplesToIgnore + 1):(end - rightSamplesToIgnore)), sweeps, span, degOfPolyFit);

return;
% end of the function "calibrateSignalsSingleMask"
